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1 Introduction 

Computational electronic Structure Theory 

Computational: 

The ever growing field of Computational electronic structure theory combines theoretical physics and 
chemistry, math and computer science. It forms an intersticial between theoretical and experimental 
physics and chemistry, gaining more and more recognition as a discipline of its own. It conqueres 
the challenges that arise from the availability of faster and novel computer hardware (e.g. massively 
paralell systems, GPU, ...). 

Electronic structure: 

The electrons are the 'glue' that keeps matter together. They determine shape and color of all things 
that surround us. On the atomic level they determine the molecular and crystal structure. On the 
level of transport properties conductivity, heat, capacity, phase diagrams and reactivity are due to 
electronic interactions. 

We distinguish between ground state properties, e.g. atomic structure, forces, reaction barriers and 
exited states which give rise to spectroscopic properties and the response of the system to external 
pertubations. 

Theory: 

It is extremly non-trivial to find a theory from firts principle, a theory that does not need additional 
parameters, to describe the electronic problem. The field is under constant and active development, 
challenging current research topics include the incorporation of finite temperatures, bond breaking 
and making, high T c superconductivity. Of special interest are strong correlations because strong 
electron-electron interactions are not described by an effective medium. 

The overall goal is to have a theory that accurately and reliably predicts material properties from 
first principles (and therefore has to be parameter free!). 
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2 The Schrodinger Equation 



For a system of electrons and nuclei the non-relativistic time-indcpcndcnt Schodingcr equation is: 

= (2.1) 

where the Hamiltonian is: 

N N M N M 

i=l i=l j>» 1 i=l a=l 1 1 a| 



M MM 
a— 1 a— 1 6>a 

position of the electrons 
position of the nuclei 
charge of the nuclei 



m e = e = h = = 1 

47re 



Ra 
Z a 

and we have used atomic units: 



The energy is given in Hartree. 

[Ha] = 27.211 eV 

It can be solved for some simple systems: 

• hydrogen atom (Laguerre polynomials) 

• harmonic oscillator (Hermitc polynoials) 

• particel in a box 

In general we have to deal with many electrons (N e i ec ~ 10 23 (!) in solids). To demonstrate the sheer 
complexity of a naive attempt of an exact solution of equation 2.2 let us consider a Silicon atom. Simply 
wanting to store the wavefunction on a grid with only 10 points in each dimensions produces a grid with 
10 3JV =10 42 points. Assuming 128-bit double precision complex numbers, each sample point requires 16 
bytes of storage, which gives approx 10 43 bytes of data. A regular blu-ray disc holds about 50 GB and we 
would therefore need ^10 32 disks - an astronomical amount. In fact, if each disk case was a centimeter 
thick the stack of disks would cover the distance from the earth to the moon ^10 21 times! 

We have to find approximations! 

Different frameworks lead to different approaches that lead themselves to different approximations that 
work for different systems amd circumstances. 
In general: 

1. wavefunction (JS) based: 

• Quantum chemistry (finite systems): Hartree Fock pertubation theory (e.g. MP2), coupled 
cluster 

• Quantum Monte Carlo (finite, periodic): variational MC, diffusion MC 
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2 The Schrddinger Equation 



2. Density based (n(r)) 

• Density-funtional theory (finite, periodic) 

- many exchange-correlation funtionals 

- exact limit not easily attainable 

- fast approximations 

3. Green's function based (G (r, r', u>)) 

• many-body pertubation theory (finite, periodic): GW, T-Matrix, BSE, FLEX, ... 

• DFT+Dynamic Mean Field Theory (DMFT) (periodic, strongly correlated) 

4. Density matrix (n(r, r') = G (r, r', t = 0)) , not covered in this lecture 

• density matrix functional theory 

• few xc approximations, little explored 



2.1 Born-Oppenheimer Approximation 

Since the nuclei are much havier than electrons they move much slower. Therefore the potential feld by 
the electrons, caused by the nuclei, can be approximated by the field of fixed nuclei. The kinetic energy 
of the nuclei is thereby neclected and the repulsion between nuclei becomes a constant. The remaining 
Hamiltonian is H e i ec (2.2), which describes the motion of N electrons in the field of M point charges. 
The electronic Hamiltonia can be solved by: 



Helec^elec ^elec^ elec 
^elec = * eiec ({r i },{R a }) 

{r^} explicitlly 
{R a } parametrically 

Additionally the total energy depends as well parametrically on positions of the nuclei. 



with the electronic wave function: 



which depends on: 



Eelec — E e i ec ({R a }) 



The total energy for fixed nuclei is: 



E, 



elec 



M M 

EE 

a=l b>a 



Z a Zb 



(2.4) 
(2.5) 



(2.6) 



(2.7) 



This is the electronic problrem we are interested in! 

Now we can reverse the argument and solve the problem for nuclei moving in an effective potential of 
the electrons. Since the electrons move much faster we can replace the electronic coordinates by average 



3 



2 The Schrddinger Equation 



values, which are averaged over electronic wavefunctions 

M , M M 



v— v 2 + vy ZaZb 

^2M a a ^^^\R a -R 

IN N M . 

(-£^ 2 +££t^-££ 



AT M 



M M 



21 t=l a =l 1 1 a 



E-nucl — 



a—1 a—1 b>a 

M 

£ ? r r V 2 + J B tot ({R a }) 



(2.8) 



The total energy of the electrones provides a potential energy surface for the motion of the nuclei. 
Forces can be calculated by the derivative of the total energy with respect to the coordinates of the nuclei. 



Etot ({Ra}) 

A 




{Ra} 



Figure 2.1: Potential energy surface 
dE tot 



F 



Hnucl^nucl = E^ nuc i = E^ nuc i ({R a }) 

^nuci ({Ra}) describes the vibrations, rotations and translations of a system. 

* ({rj , {R a }) = V elec ({r. ( } , {R Q }) V nucl ({R a }) 



(2.9) 
(2.10) 

(2.11) 



is the full wave function of the system. 

The adiabatic approximation fails for system where non-adiabatic effects, e.g. electron-phonon coupling, 
are import. The describtion of superconductivity is not possible in this context. The Born-Oppenheimer 
approximation also fails in cases where more than one energy landscape is present (molecule approaching 
surface). 



The central goal of this lecture is to solve the electronic problem 

Helec^elec = E e i ec §> e ; ec 

which is given in (2.2). 
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2 The Schrddinger Equation 

2.2 Antisymmetry or Pauli Exclusion Principle 

Electrons are Fermions and therefore carry a spin. The Schrodinger equation is not explicitly spin depent. 
We therefore have to introduce spin functions in spin variable u. 

t : a(a) 
I : 0{<t) 

These have to be complete and orthorgonal (for convenience). 

J daa* (a) a(a) = J dcr/3* (a) /3 (a) = 1 

J daa* (cr) (3(<t)= J dcr/3* (a) a (a) = (2.12) 

The electrons arc described by a collective coordinate: 

x{r,<7}, (2.13) 

the wave function is: 

* = ({xi}). (2.14) 

Antisymmetry: 

A many-electron wave function must be antisymmetric with respect to the interchange of the 
coordinate x of any two electrons. 

*(xi,...,x i ,...,x J -,...,xj V ) = *(xi,...,x J -,...,x i ,...,xj V ) (2.15) 

(Generalization of the Pauli exclusion principle) 
The many electron wavefuntion must satisfy: 

1. The Schrodinger equation 

2. Be antisymmetric. 
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3 Wave function based approaches 



Let us imagine we only had one electron or we could postulate an effective electron-electron interaction 
that describes the effect of all other electrons in an average way. 




JV , N M 



2 f-> ^ |r 4 - R a 

z— 1 2—1 a— 1 1 



ft. (i) (one electron) 
AT M 



Ti — I\ 

1=1 j>i 



couples all electrons^ makes life difficult 



3.1 Hartree theory 

In quantum mechanics we have the wave-particle duality; namely the notion that the electron's position 
in space is given by a probability distribuntion. 

n(r) = |^(r)| 2 (3.1) 

we call tp(r) the spatial orbital or wave function for the electron 

|^(r)| 2 dr 3 (3.2) 

is then the probability of finding the electron in the small volume element dr. Combinded with its spin 
we can define spin orbitals: 

or (3.3) 
^(r)/?(<r) 

To make progress we next consider a system of N non-interacting electrons having a Hamiltonian: 

N 

H NI = J2Hi) (3-4) 

i=l 

h(i) will general include the kinetic and potential energy of electron i, but may also contain an effective 
electron-electron interaction. 
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3 Wave function based approaches 



Each operator h(i) will have a set of eigenfunctions that we take to be a set of spin orbitals. 

h(i)tpj (xj) = ejipj (xj) (3.5) 

which are orthorgonal for convenience: 

J dxtfi (x) ipj (x) = Sij (3.6) 

What is the corresponding cigenfunction of H then? 
Answer: 

a product of the spin orbitals 

V HP (xi, . . . ,xjv) = fi (xi) ipj (x 2 ) ...<pk (xjv) (3.7) 

The eigenvalue of 

jjNI^NI = m NI (3 g) 

is simply the sum of eigenvalues, because h{i) only acts on the orbital with corresponding coordinate x^: 

E = ei + ej + . . . + e k (3.9) 

The eigenfunction ^ HP is alos called a Hartree product and is an example for an uncorrected wavcfunction, 
because the probability of finding electron 1 in rfxi, electron 2 in dx.2, etc., is simply the product of 
probabilities: 

\ipi (xi)| 2 \dxwj (x 2 )| 2 dx 2 ... \(fi k (xat)| 2 dx N (3.10) 

or, if we view this differently we have factorized the wave function; a concept that is recurring throughout 
physics. 

So if we forget for now that ty HP does not obey the antisymmetry principle we could make the ansatz: 

H* HP = E^ HP 

and ask, which single particle orbitals minimize the energy under the constraint that the orbitals are 
normalized. 

=> Lagrangian: L[V] = E [*] ~J2 €i dx< ^* ( x ) W ( x ) 

i J 

The energy is given by 

because ty HP is normalized. 
We already know that: 

<*™>(i)l** P > = (<Pi\h(i)\<Pi) 

For the Coulomb potential we get: 



HP 



1 



.HP 



[ dXj . . . dXjV<£* (xi) ...<pl (xjv) i — - — Tfn (xi) ...<p k (xjv) 
J Tj r,- 



/ / ' " ' ' ' V. _ , 

3 I 

dXidx . <Plte)<P* m ( X j)y» ( X »)^ ( X j) 



All other combinations of ipk integrate out 
I and m are the occupied orbitals 
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3 Wave function based approaches 



HP 



N N ^ 
^— < ^— < | r j.' 



HP 



Now our Lagrangian reads: 

N N N 

L [$] = ^ / dx^i (xj) ft (xj) (fi (xj) = y^y^ / dxdx' 



gf |^w>iwyww (3 . u) 



yVi ( X ) ¥>j( X ')^ W^(X') 



r - r' 



doing the varition with respect to <p* and demanding: 

£=0 V, 

gives: 

= ft (xO ^ (Xi) + ^ / dx ^ X V , (x) - ^ (x) 

j>i ■'II i 



/^R*W-/^W 

S v ' 

V H (r) ipi (x) . 

[ft (r) + V H (r)] ^ (r) - J dr'^L^ (r) = (r) (3.12) 

are the Hartree equations. 
A few points to note: 

• Vh, the Hartree potential depends on the electron density: Vh = Vff(n(r),r) and the density 
depends on the orbitals. Therefore the Hartree equations have to be solved self-consostently. 

• The sum in Vh runs over all electrons and therefore contains the interaction of an electron with 
itself. This is removed exactly by the 2nd, orbital dependent term. 

• ft (r) and Vh (r) are the same for all i electrons but the self-interaction correction introduces an 
orbital specific dependence. 

• The Hartree potential is the classic electrostatic potential of a charge distribution, it is positive and 
therefore repulsive, e.g. for the hydrogen Is function: 

i>u (r) = ^e _r 

we get 

Mr) = i(l-(l + r)c- 2 ''), 

which looks like: 
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3 Wave function based approaches 



V H (r) * 




1 

r 

Figure 3.2: Hartree potential 

It is repulsive and delocalizes the electrons. 
• The Hartree potential largely counteracts the external potential: 

\" \ " Z a ST f HP n " ^ 

2^| ri -B,| If} dK |r,-R| 

V es , t (r 4 ) 

3.2 Hartree and Hartree Fock 

ty HP is an eigenfunction of the many-electron Schrodinger equation, but it does not satisfy the antisym- 
metry requirement. 
As is easily seen: 

V HP (xi, ...,Xi,...,Xj,.. . ,xjv) = (xi, . ..,Xj,.. . ,Xj, . . .,xjy) (3.13) 
Considering the two-electron case - the hydrogen molecule: 




Figure 3.3: Two electron case 



*12 ( x i> x 2) = <pi (xi)<Pj (x 2 ) 
*?i P (xi,x 2 ) = ipi (x 2 ) </Jj (xi) 
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3 Wave function based approaches 



The two wavefunctions clearly distinguish the electrons. 
Indistinguishable and antisymmetric: 

1 



* (xi,x 2 ) = -j=[fi (xiVj (x 2 ) - <pi (x 2 )v?j (xi)] 

Via d - c b 



(3.14) 



This is mathematically equivalent to a determinat: 

A=( a b ) detA = 



a b 
c d 



= ad — cb 



*(xi,x 2 ) = 



1 

71 



tpi (xi) ipj (xi) 
<Pi ( x 2) <Pj (x 2 ) 



t electrons 



orbitals 



or in general 



$ SD (x 1; ...,x w ) 



1 



ipi (Xl) Ifj (xi) ... ip k (xi) 
(x 2 ) <fj (x 2 ) ... ^ fc (x 2 ) 



(3.15) 



<Pi (Xjv) ¥>j (Xjv) • • ■ ¥>fc (Xjv) 

and is commonly called a Slater determinant. 

Next we will approximate the many-electron wave function by a single slater determinat and use the 

variational principle to determine the spin orbitals that minimize the grounnd state energy. 

This gives the Hartree-Fock method/approach. 

First we need the energy for a Slater determinant (SD): 



E 



(^ SD \H\^ SD ) 



As we have seen H e i ec has two distinct parts: 



(* SD \H\* SD ) 



~>& SD is normalized 



N N 



(3.16) 



i i j>i 

We therefore need to evaluate two different contributions: 



(^ bu \h(i)\^ blJ ) 



SD\ 



1 



\* SD ) 



For simplicity this will be illustrated for the two-electron example: 

J dxidx 2 -^=[<# (xi)(pj (x 2 ) - tfi (x 2 )^j (xi)]/i(l)[< / j i (xi)^ (x 2 ) - tfi (x 2 ) <pj (xi)] 
= y dxidx 2 -^=[vi (xi) ipj (x 2 ) ft(l)v?i (xi) ^ (x 2 ) - <pi (xi) <pj (x 2 ) h(l)ipi (x 2 ) ^ (x x ) 



+ <fii (X 2 ) </?j (Xl) /l(l)<£j (X 2 ) ifj (Xl) - ^ (X 2 ) (Xl) /l(l)<£j (Xl) Ifij (x 2 ) 



After integrating over x 2 : 



/ dXl 7f ^ (xi)/l(l)Vi (Xl) - Ifj (xi)/i(l)^- (xi)] 
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3 Wave function based approaches 



h{2) yields an identical result. 

(V SD \J2 h^ SD ) = (xi) IMObi (xi)> 

This is the same result we obtained in the derivation of the Hartree approximation. The general proof can 
be found in Szabo and Ostlund (and might be given as an appendix). 
For the Coulomb operator we again consider two electrons first: 



J j dxidx 2 i [ip* (xi) 95* (x 2 ) - if* (x 2 ) ¥?* (xi)] | r . ^ 1 

* [<Pi ( x l) ( x 2) - <Pi (x 2 ) V'j ( x l)] 

= (( dxidx 2 i 1 1 

*[ a) tfi (xi)^j (x 2 )</?i (xi)<£j (x 2 ) 12 12 

6) (x 2 ) ^ (xi) ^ (x 2 ) <fj (xi) 21 21 

c) -</?i (xi)v?j (x 2 )<^i (x 2 )^- (xi) 12 21 

d) -<fi (x 2 ) <fij (xi) <^ (xi) <pj (x 2 )] 21 12 

4* 4' *r 4" 

• down the columns (!) always the same orbitals 



• across (— >•) filled by different combinations 

However, this is not the most convenient way of writing it. So let's reorder each term so that we have a 
different integration order. 





Xl 


x 2 


Xl 


x 2 


a) 


i 


j 


i 


j 


b) 


j 


i 


j 


i 


c) 


i 


j 


j 


i 


d) 


j 


i 


i 


j 



We obtain the same pattern, but now in the orbital indicies. 
If we define: 

(ij\kl) = J dxdxVi (x) ifj (x') - ^ — ip k (x) ipi (x') 

then the matrix element for the Coulomb operater becomes: 

1 



SD 



SD 



or in general with: 



N N 



V lr- - r I 

i j>i 1 1 3 < 
N N 

^sd |yee| ySD^ = - ^ ^ (mn\mn) - (mn\nm) 



m n=tm 



(1) 



(2) 



(3.17) 
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3 Wave function based approaches 



(1) As we will see, this will again give rise to the Hartree potential 

(2) This term is new and gives exchange because indices are exchanged 
=>■ The energy of a single Slater determinant: 

N N 

2 



with 



Ehf = (V SD \H elec \ * SD ) = ^2(n\h\n) + - ^(mn|mn) - (mn\nm) 

n nm 

(n\h\n) = J dxtp* n (x) h(x)(p n (x) 



Now we make one important observation: 



Altough we have a many electron wave-function and a many electron hamiltonian we can write the 
energy of a Slater determinat in terms of 1 and 2 electron integrals only. 



=>■ This of course follows from the fact that the Coulomb potential is a two-electron operator! 

=> As we will see later: This is one of the key foundations of quantum chemistry 
To find the minimum of the Hartree-Fock energy we again define a Lagrangian: 

N 



L [0> SD ] = E [V SD ] - £ e m [ dx^4 (x) Vm (x) 

N . N . 

= E / dx ^™ ( x ) M x )<£m (x) - ^ e m / dx<^„ (x) ip m (> 



N N 

+EE / dxdx ' 



<Pm ( X ) <Pn ( X ') <Pm (x) (x') ip* m (x) (f* (x') <p n (x) ifi m (x') 



r - r' 



r - r' 



(3.18) 



• kp* m shows up twice in (3.18) =4> cancels factor of 2 
and in analogy: 



SL 

6<fi* 



= Vi 



A" 



h(x)<p m (x)+^2 J dx ' 



'fin ( X ') ¥>m (x) (/?„ (x') ip* n (x') (x) lfi m (x') 



r - r' 



r - r' 



(x) / dx' 



(x')Vm(x) Pm(x') 



= e m (^ m (x) 



[fc(x) + y H (x)] ^ m (x) - E / dx' ^" |^" (X) (x') = e m y m (x) (3.19) 



Fock Operator / 
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3 Wave function based approaches 



To rewrite equation (3.19) in a more convenient form we apply from the left J2 m I ^ X( Pm ( x ) : 

£ e m = W dx^;„ (x) [Mx) + V H (x)] Vm (x) - J2 I dxdx ^- (x/) ^ /)y " (x) yTO (x') 

m m nm ' ' 

TV N 

— 5^( n l^l n ) + y^(mn|mn) — {mn\nm) 

n nm 

l N 

E HF = — y^(mn|mn) — (mn\nm) 

1 * 

£ flF = ^ e m — - y^^(mn\mn) — (mn\nm) (3.20) 



2 ■ 

nm 



(*) 



(*) cancels double counting of Coulomb repulsion and exchange in the sum of eigenvalues (which counts 
interaction i o j once for i and a 2nd time for j) 



we define the density matrix: 



.(x,x') = 5>;(x)^„(x') (3.21) 



(note: the diagonal of the the density matrix is the density) 
We also define the non-local potential: 

^(x,x') = -^^ (3.22) 

=> [Mx) + V H (x)] <p m (x) + J dx'S HF (x, x') Vm (x') - e m <p m (x) (3.23) 
Often this equation is also written like this: 

h(x)<p m (x) + J dx'V HF (x, x') tp m (x') = e m tp m (x) 

I I 
independent of {ipi (x)} dependent of {ipi (x)} 

f = h + V HF is the Fock operator (3.24) 

• The Fock-operator is an effective one electron operator 

• It replaces the many-body Schrodingcr equation by a set of one electron equations, where each 
electron moves in an effective field, also called mean-field 

• The mean-field is an important concept because it allows us to seperate the many-electron problem 
into one-electron problems: 
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3 Wave function based approaches 




Figure 3.4: Effective Potential 



• The wave-function of electron m does no longer explicitly depend on that of all other electrons only 
implicitly through Vrf 

• In fact, it is actually possible to speak of one-electron wavefunctions or orbitals, which is not really 
possible for the full many-body wave- function (xi,X2,X3, . . . , xjy). 

• The Fock operator has two parts 



/ = 



HF 



{ifi (x)} -independent 
>It has to be solved self-consistently. 



V 
X 

{ifi (x)} -dependent 



V HF is called the self-consistent-field (scf) in which the electrons move. 

1. Start with an initial guess for {(f t} 

2. Calculate n (x), n (x, x') and V HF (x) 

3. Solve the HF equations (3.19) for a new set of {fi] 

4. Check for convergence (is n (x), n (x, x') or the total energy the same, within a given tolerance, 
as in the previous iteration. 

- If not: repeat calculaion (go to 1.) 



- If yes: finished 

M r \ l \ [ A ' \ " fn ( x ') <Pn ( x ') / n 

• V H (x) ip m (x) = J dx | r -r^| Vm ^ 

• J dx S (x, x ) ip m (x) = J dx | r ~ r ,| ¥n (x) 

<fim have been exchanged arises from antisymmetry, S is called the exchange operator 

Hartree Theory Vh' Coulomb 
Hartree-Fock Theory: Coulomb and echange 

Befor we proceed to an expression of the Hartree-Fock equation that one would solve in the computer let's 
have a different look at the Coulomb and exchange terms. 
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3 Wave function based approaches 



3.2.1 Closed shell Hartree Fock 

Let's consider a system with an even number of electrons and this scenario: 



N/2 



- # 2 

Figure 3.5: Closed Shell 

The spatial orbitals are restricted to be the same for spin j" and J,. 
==> we can pair up orbitals: 

ipi (x) = -01 (r) a (u) = -01 (x) 
<p 2 (x) = 0i (r) (3 (oj) = (x) 



=>• The Slater Determinat of tp±, (f2 
and the energy: 



yjv becomes a Slater determinant of the functions ipi , ipx, 02, 02, 
E 1 = ^~^(n|ft|n.) + — (mn\mn) — (mn\nm) 

n n n^m 

N N N 

^~^{n\h\n) H — ^^^^(mn|mn) — (mra|nm) 



0JV/2, 



JV 



iV iV 



N/2 

Y,(n\h\i 

n 



N/2 



N N 



mnlmn) — (mnlnm) 



drda0; (r) a* (a) ft (r) 0„ (r) a (a) - / drda0; (r) /3* (a) ft (r) 0„ (r) /3 (a) 
dr0* (r) ft (r) 0„ (r) = (n\h\n) = h„, 
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3 Wave function based approaches 



N/2 ^ N N 

= 2 y~^(n|/t|n) + — y^(mn|mn) — (mn\nm) 



n m 

N/2 N/2 N/2 N/2 

^2 ^2 PmVn = ^2(¥>m + <Pm) ^(fn + <fn) 
m n m n 

N/2 N/2 



m n 

N/2 N/2 N/2 



— 2 'y^(n\h\n) + — ^^(mn|mn) — (mn\nm) 



2 

n m 

+ (mn\mn) — (mn\nm) 
+ (fhn\mri) — (mn\nm) 
+ (mn\mn) — (fhn\nfh) 
aa'au' aa'aa' 

dada' . . . 

For these integrals to be non-zero the 1st and 3rd position and the 2nd and 4th position have to both have 
a bar or non. 



>with (mn\kl) ^/drdr ^"^^^'^ 



N/2 N/2 

E = 2 y^(n\h\n) + V2 (mn\mn) - (mn\nm) (3.25) 

n nm ^ 



r drdr , fc(r)^(r')i(r)i(r') = f ^ (r)| 2 (r')| 2 
J |r-r'| J |r-r'| 



2 



The Coulomb integral J mn is the classical Coulomb repulsion between two charge distributions \ip m (r')| 
and \tj} n (r')| 2 

J |r-r'| 
The exchange integral K mn has no classical analog. 

3.2.2 Hartree Fock in real-space 

The Hartree Fock Self Consistent Field Equation in real-space follows either from the expression for E HF 
by minimizing the Lagrangian or by integrating out the spin in the spin-dependent HF-equations: 

N/2 

f (r) V™ (r) = [h (r) + V H (r)] ^ m (r) + J dr' £ ^Mj^l ^ (0 - ^ m (r) (3.28) 

V v ' 

/" dr'£(r,r')Vm(r') 

The kinetic energy can be written as high order finite difference expansion on a grid (see fig. (3.6)). 

M 

V 2 ^ = ]T C n 4>{x t + nh,y^z k )+0(h 2N+2 ) (3.29) 



n=-M 



For further information see CD. Smith - Numerical solutions of Partial Differential Equations, Finite 
Difference Methods, Oxford University Pree, new York, 1978. 
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h 



A 



y 





X 

Figure 3.6: Integration grid 

With this expansion for the kinetic energy the Hartree fock equations could be solved on a grid with 
suitable quadratures for the integrals. 
However, the pproblem are the core states! 

Take e.g. a Is state (fig. (3.7)): It is tightly bound to the nucleus and many grid points are required to 
resolve the wavefunctions near the nucleus. 




nucleus / core 

Figure 3.7: Is state radial function. 
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nucleus / core 



Figure 3.8: 3s state radial function. 



Even a 3s state (fig. (3.8)) has rapid oscillations in the core region to satisfy orthogonality. 

3.2.3 Hartree Fock in a basis 

Solution to this problem: 

• Introduce a basis 

• build the rapid ocillations in the core region into the basis e.g. basis functions could be orbitals of 
free atoms 

Without specifying the type of basis functions yet we introduce a set of K-many basis functions. 



The molecular orbitals can then be expanded in this basis: 

K 

i>i 0) = ( r 

• If {4>^} was complete this would be an exact expansion 

• In general {(f)^} is not complete 

Inserting the expansion into the Hartree Fock equation (3.19). 



(r), M = !...#} 



(3.30) 



with 



/ = F (r, r') =h(r)5(r- r') + V HF (r, r') 



we obtain 




K 



K 



(3.32) 
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Acting with J dr</>* (r) from the left yields: 



K . K ,. 

YtCn J drdr'^F(r,r')^ (r , ) = e i ^C Mi J dr# (r) ^ (r) 



(3.33) 



Fffl S U fj 

Overlap matrix S vfl : Basis functions are in general not orthonormal 



K 



K 



(3.34) 



• We have transformed the Hartree Fock equations for single particle states to a Hartrce Fock equation 
for the expansion coefficients 

This can be written in matrix form. 

FC = SCe (3.35) 

Where C is a K x K square matrix of the expansion coefficents and e a diagonal matrix of singe particle 
energies. 



e = 



\ 



£2 



(3.36) 



V e K / 

Note that equation (3.35) produces only as many eigenvalues as there are basis functions, i.e. K-many. 
==>- K-many basis functions can only expand K-many single particle states. 
Equation (3.35) is the Roothan equation, which is central to quatum chemistry. 
We have moved form: 

Helec^ elec — Eelec^& elec 

to the Roothan equation: 

FC = SCe 

in two steps: 

• Single dctcrminat approximation (conceptual or physical approximation) 

• Introduction of a basis (expansion of wavefunctions) (practical or computational approximation) 

To solve Roothaan's equations all that is left to do is to derive equations for F in the basis and for a given 
basis term this into equations that can be implemented. 

F(r,r') =h (r) + v H (r) + £ (r, r') 

i 



(3.37) 



The density: 



N/2 N/2 

.(r)=2j>i (r)| 2 = ][>*(r)^(r) 

i i 
N/2 

= E E 2<7 - C ^ €(r)<Mr) 



P M , y - density matrix 
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3 Wave function based approaches 

The Hatree potential: 

V* = J dr# (r) v H (r) ^ (r) = J drdr'^ (r) (r) 
/ r — r' 

nm 1 1 

nm 

The exchange operator: 

= | drdr'C (r)£(r,r')«^ (r') 

= /drd r ^(r)E^§P^(r') 

= / drdr'C (r) £ £ ^™ ^-^ ^ ^ 



~P n m{vn\mfi) 



J? 



h vfl + E {vn\u.m) - ^(vn\mn) 



nm 

independent on P 



dependent on P 

The Roothaan equations are non-linear. 

F = F(P) = F(C) 
F(C)C = SCe 

Iterative or self-consistent solution of the Roothaan equation: 

1. Choose a basis 

2. Calculate the integrals £ MI/ , h uv , (un\fj,m) 

3. Initialize P nm 

4. Compute F uv 



5. Solve the Roothaan equations New coefheients C v i => new P, 
Repeat steps 4 and 5 until P nm does not change anymore. 



nm 



(3.38) 



Scaling: 

For an atom centered basis it is easy to see that our basis {(j)^} grows liniear with the number of atoms N. 
=> TV 4 many integrals {vn\^m) are required 

=>■ Computation of F^ v requires iV 4 operations if no further tricks and reductions are employed 
=>■ The formal scaling of Hartree Fock is iV 4 
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Figure 3.9: Gaussian 



One solution, a fast basis: 
Gauss orbital (or Gaussians) 

"0n 7 (r) = A» x n v 7» , z" 1 y n " z"* e~ 7T " 2 Kartesian Gaussians (fig. 3.9) (3.39) 

with n x + n y + n z — L we obtain the familiar functions from spherical symmetry (tab. (3.1)). 
Atom centered Gaussisans: 

</>an7 (r) = (r — R a ) with R a the position of atom a (3.40) 

Important properties: 

• Fourier transformation of a Gaussian is also a Gaussian 




(3.41) 



benefitial for solid state calculations 



L 


# 


orbital type 


n x 


n y 


n z 


n7 (,(f) 3/ V- 2 ) 





1 


s 











i 


1 


2 


Px 


1 












3 


Py 





1 









4 


Pz 








1 




2 


5 


d xz 


1 





1 


A'yxz 




6 


dyz 





1 


1 





d xy 110 A-yxy 

d 3z2 _ r , (3z 2 - r 2 ) 



9 d x 2_ y 2 27 (x 2 - y 



10 



■2\ 



7TI 



Table 3.1: Parameters for cartesian Gaussians up to n x + n y + n z = L = 2. The last three rows for L 
follow from the linear combination of x 2 , y 2 , z 2 -orbitals. 
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The product of two Gaussians yields another at a different position: 



c - Q (r-R )V 



-Rb)" = e -5-D^ e -«(r-X) 



with the center of mass X = aR °+^ Rb and D = R a - R fc as well as, S = -^Ar. 

a+p u u 5 a+p 

For the cartesian functions a similar relation can be derived (fig (3.10)): 



(3.42) 



f3(x-R b y 




e -5D' e -£(x-X) 



Figure 3.10: Product of two s-type Gaussians 



i+j 



( X -R a y(x-R b y = EE 



n=0 I 



" 1 ' D l+ i- n {x-X) 



n-l) {/)' ^ in - /i" ' " 



T ijn = ^ i+ J- n e- 5D2 [] Mi, 



(3.43) 
(3.44) 



i+j 



ia (r) (pbjfl (r) = ^ Tijn</>nX£ (r) 



(3.45) 
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X 




Figure 3.11: Product of Gaussian orbitals on two different atoms 

• Products of Gauss orbitals on two different atoms can be expanded in Gauss orbitals centered on 
their center of mass 

• The transformation tensor is given analytically, to compute k center integrals 

jtl (r)ft (r')^ (r)^ m (r') 



{vn\[im) — J drdr' 



(3.46) 



we have to expand the two products for r and r' into Gaussians around their respective centers of 
mass. 

Then the problem is reduced to calculating the Coulomb integral of two Gaussians. 





Figure 3.12: k-center integrals 

Let us consider two s-type Gaussians at the centers R 9 and R p . We then have to solve the following 
integral (the derivation is quite tedious and can be done by inserting the Fourier transformations of 
the Gaussians and p): 



drdr'- 



o\r— R p | 



Pi e 



Flo 



27rf 



— ; U2 F ° ( — i — l-^p _ ^-9 1 2 

pqip + l)' \P + Q 



(3.47) 
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with 

F (t)=t^f\y C ^= 1 -^) 1/ \r i (t^) (3.48) 

erf is the error function, avaiable in any/most programming languages. 
The Coulomb integral of two Gaussians is known analytically. 

Now all that is left to do is to determine the Cpolomb integral of two arbitrary cartesian Gaussians. We 
do this by observing that the derivative if an s-type Gaussian produces higher order Gaussians. 

for one dimension <f> n (r — R) = (r — R) n e~ a ^ r ~ R ^ 

t 

note: normalization is different 

^° ( ; ~ R) =2a(r- R) = 2a ^ (r R) 

aR 

d20 °J;r R) = -2ac-^ 2 + 2a(r- R) d *> ^ ~ R) 
AR AR 

= 2a l-fa (r-R) + 2a (r - R) 2 e^"^ 2 
= 2a [2a<j) 2 (r - R) - <fo (r - R)} 

So we make the ansatz: 

^^^ ^Ai^-fe-^ (3.49) 

j<i 

Which gives us the recursion relation: 

Da = -2aA-i, j-i + U + l)A-i, j+i ( 3 - 5 °) 
with D o = 1 and D io =0 V i,j<0Aj>i 

=> Since derivations in x, y,z, are seperable: 

(-["a: jn B (\n z mi<nj 
V£,0 aO a(r) = -^---^--^-^(r) = fnm^amc,(r) (3.51) 

x v z in 
The diagonal elements of D nm arc non-zero and therefore D nm is invertable and we can write: 

amQ (r) = J2 ^„mV^ao a (r) (3.52) 



=> In other words, the other cartesian Gauss functions follow from the s-functions through a recursion 
relation involving the derivatives. Since the Coulomb integral of two s-functions is known analytically (see 
eq. (3.47)), as well as their derivatives Vj^ a and Vj^ b (see eq.(3.51)), the Coulomb integral of 2 arbitray 
Gaussians can be computed fast via a recursion relation. 
Recipe for computing (vn\fj,m) (also see fig. (3.12)): 

1. Expand products 4> v <pfj, and 4> n 4> m in Gauss orbitals around their centers of mass 

2. Evaluate recursion relations for the Coulomb integrals on the centers R p and R 9 
Advantages: 

• No 6-dimensional real-space integral J drdr' has to be performed => fast 

• Efficient algorithms and screening prosedures have been developed for carrying out the expansions 
and recursions fast, even on the fly evaluation of 4-center integrals is possible =>■ on the fly 
saves memory! 
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01s 



Figure 3.13: Comparison between Gaussian Is function and actual Is function 

Remark: 

• The procedure outlined here is for illustrative purposes, but by no means the most efficient. 
Disentvantages: 

• Gaussians miss the cusp for Is functions (certainly not ideal for core states) 

• Difficult to design basis sets (exponents, number of functions) that allow for a systematic accuracy 
analysis 

3.2.4 Performance of Hartree-Fock 



Table 3.2: Hartree-Fock total energies for selected atoms and simple molecules. 



Eft [Ha] 


H 


He 


Be 


Mg 


H 2 


LiH 


HF* 


-0.5 


-2.8610 


-14.5697 


-199.6900 


-1.1371 


-7.9872 


cit 


-0.5 


-2.9362 


-14.6674 


-200.0530 


-1.1660 


-8.0400 


deviation (CI-HF) 


0.0 


-0.0427 


-0.0977 


-0.3630 


-0.0289 


-0.0528 


in % 





1.5 


0.7 


0.2 


2.5 


0.7 



*FHI-aims calculations with cc-pV5Z Gaussian basis set 
tfrom Chakravorty et. al. Phys. Rev. A 47, 3649 (1993) 
L. Wolniewicz J. Chem. Phys. 99, 1851 (1993) 
X. Li and J. Paldus, J. Chem. Phys. 118, 2470 (2003) 



• Hartree-Fock underestimates the total energy compared to escentially exact Configuration Interac- 
tion calculations 

• The underestimation is only a small persentage of the correlation enrgy 

• However, this can be decisive for energy differences 
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-0.7 



ao 
a 

oj ■ 

3 
o 
H 




o-o Restricted HF 
— Exact 



2 3 
Bond length (A) 

Figure 3.14: Total energy of H2 in restricted closed-shell Hartree-Fock as a function of the interatomic 
distance. 



Figure (3.14) shows the restricted Hartree Fock energy for H2 as a function of the interatomic distance. 
The curve shows 3 distinct regimes: 

- Repulsion at short distances 

- A minimum at intermediate distances (to see if it actually binds one has to look at the binding energy) 

- Levelling out at dissociation 



o 



(a) Wavcfunction of the bonding state of H2 at equilib- 
rium in restricted closed-shell Hartree-Fock. The two 
spheres mark the atoms, the ellipsoid is an isosurface 
plot of the orbital and the shading corresponds to the 
wavefunctions on a 2D cut through the molecule. 




(b) Isosurface of the wavcfunction for the anti-bonding 
state of H2 at equilibrium in restricted closed-shell 
Hartree-Fock. The lobe on the left corresponds to nega- 
tive values and the lobe on the right to positive ones. 



Figure 3.15: The wavefunction for the bonding and anti-bonding state of H2 at equilibrium in restricted 
closed-shell Hartree-Fock. 



Figures (3.15, 3.16) show the two lowest Hartree Fock single particle orbitals at the equilibrium bond 
distance and close to dissociation at 6A. 

• The bonding state is occupied by 2 electrons 
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(a) Wavcfunction of the bonding state of H2 at a binding 
distance of 6 A in restricted closed-shell Hartree-Fock. 
The two spheres that mark the atoms are hidden under 
the isosurface plot of the orbital. The shading corre- 
sponds to the wavefunctions on a 2D cut through the 
molecule. 

Figure 3.16: The wavefunction for the bonding and 
in restricted closed-shell Hartree-Fock. 




(b) Isosurface of the wavcfunction for the anti-bonding 
state of H2 at a binding distance of 6 A in restricted 
closed-shell Hartree-Fock. The sphere on the left corre- 
sponds to negative and the one on the right to positive 
values. 



i-bonding state of H2 at a binding distance of 6 A 



The antibonding state is empty 

The bonding state derives from a positive superposition of the Is states of the two hydrogen atoms 
and the anitbonding state from a negative superposition (This is most obvious at dissociation (6A)). 



0.3 




0-0 Restricted HF 
— Exact 



2 3 o 

Bond length (A) 



Figure 3.17: Binding energy of H2 in restricted closed-shell Hartree-Fock as a function of the interatomic 
distance. The exact curve is taken from L. Wolniewicz, J. Chem. Phys. 99, 1851 (1993) 

Figure (3.17) shows the binding energy 



Ebind = E (H 2 ) — 2E (H) 



(3.53) 



compared to exact calculations. 
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-20 Q 1 1 1 1 1 1 1 1 1 1 1 1 1 

3.0 3.5 4.0 4.5 ^ 5.0 5.5 6.0 

Bond length (A) 

Figure 3.18: Binding energy of Ar2 in restricted closed-shell Hartree-Fock as a function of the interatomic 
distance. The reference curve bt K. T. Tang and J. P. Toennies (J. Chem. Phys. 118, 4976 
(2003)) is based on theoretical modeling of experimental data. 

• The equilibrium distance in Hartree Fock is good 

• The binding energy is underestimated 

• The dissociation limit is incorrect 

Figure (3.18) shows the binding energy curve of a van der Waals bonded dimer. 

• Ar atom is closed shell - no covalent bond 

• Binding energy is much lower (mcV) 

• Hartree-Fock gives no binding in this case! 

• But correct dissociation limit 

Theoretical descrition of bond breaking and bond making is a great challenge! 
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Table 3.3: Individual energies (in Ha) for reactants and the transition state of the reaction 
OH + H2 — > H + H2O. The deviation is measured with respect to accurate 
Coupled Cluster (CCSD(T) CBS cc-pV(Q,5)Z) calculations. The experimental 
value is taken from Y. Zhao and D. G. Truhlar, J. Phys. Chem. A 110, 13126 
(2006). 



OH 


H 2 H 


H 2 


TS 


HF -75.4284825 


-1.1336606 -0.5000000 


-76.0680526 


-76.5283317 


deviation (CCSD(T)-HF) -0.2482897 


-0.0409171 0.0000000 


-0.3092634 


-0.3148024 


reaction barrier forward 


reverse forward 


reverse 




HF 0.92 eV 
Exp. 0.22 eV 


1.08 eV 21.2 kcal/mol 
0.92 eV 5.09 kcal/mol 


24.9 kcal/mol 
21.2 kcal/mol 





Hartree Fock in chemestry (see tab. (3.3): 

• Typically underbinds by about AeV per bond =>■ too large for thcrmo chemistry (energy for disso- 
ciation) 

• lkcal/mol (~ l/20eV) error in transistion state energy of a chemical reaction gives a factor of 5 
error in the reaction rate 



Table 3.4: Lattice constants (in A) of some solids in HF. Values are taken 
K. Doll and H. Stoll, Phys. Rev. B 56, 10121 (1997) 
cohesive energy (eV) LiF NaF KF LiCl NaCl KC1 



HF 


4.011 


4.636 


5.450 


5.281 


5.791 


6.548 


Exp. 


4.010 


4.609 


5.311 


5.106 


5.595 


6.248 


deviation 


0.02% 


0.6% 


2.6% 


3.4% 


3.5% 


4.8% 



Table 3.5: Cohesive energy (E(solid) — E(atonii)) of some solids in HF. Values are taken from R. Dovesi et 
al. Reviews in Computational Chemistry, Volume 21 and K. Doll and H. Stoll, Phys. Rev. B 56, 
10121 (1997) 



cohesive energy (eV) 


LiF 


LiCl 


NaF 


NaCl 


KF 


KBr 


KC1 


MgO 


Si 


Be 


HF 


6.81 


5.70 


5.95 


5.31 


5.65 


4.94 


5.54 


7.25 


6.67 


1.87 


Exp. 


8.92 


7.24 


8.00 


6.68 


7.70 


6.27 


6.75 


10.26 


9.50 


3.31 


deviation 


23.7% 


21.2% 


25.6% 


20.5% 


26.6 


21.2% 


17.9% 


29.3% 


29.7% 


43.6% 



Hartree Fock in condensed matter (see tab. (3.2.4)): 

• Not often applied 

• Lattice constants overestimated 

• cohesive energies underestimated significantly (like in chemistry) 
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3.2.5 Koopmann's Theorem and the meaning of Hartree-Fock eigenvalues 

Eigenvalue equation: 

J dx'F(xL,x')tp m (x') = e m <p m (x') | x J dx<£* (x) 
J dxdx'ipl (x) F (x, x') ^ m (x') e m S nm 



(p n (x') = e r 



- / dxdxv; (x, [» ,x) { (x - x', 1 2ii^> + £ i££laM 

L i i 

AT 

=>h nn + y~](ni\ni) - (ni\in) = e n 

i 

Now consider an occupied state (a): 

AT 

e a = h aa + ^2(ai\ai) - (ai\ia) 

i 

N 

= h aa + ^(ai|a«) — (ai\ia) 

because (aa\aa) — (aa\aa) = 

e a is the kinetic and external energy of an electron in state a plus Coulomb and exchange energy with 
all remaining TV — 1 electrons. This suggests that e a is the removal energy of electron a or the ionization 
potential (IP), a is the highest occupied state. 

IP = E(N — 1) — E(N) (3.54) 

in Hartree-Fock: 

E(N-l) = (*f r D _ 1 \H elec 9f f 1 i 1 ) 
E(N) = (* s N D \H elec * s N D ) 

SSfff 3 and ^f f D _ 1 are generally not composed of the same orbitals. If for a moment we assume that they 
are, we obtain: 



E(N) = y^(a\h\a) + ^ y^(ai\ai) - (ai\ia) 

a ai 
N 1 N 

E C (N - 1) = y^(a\h\a) + - ^ (ai\ai) - (ai\ia) 



2 



where c is the electron that has been removed. 



IP = E(N - 1) - E(N) 

1 N 1 N 

= -( c \ h \ c ) - 2 X ( ac l ac ) _ ( ac l ca ) ~ 2 ^ ~ 

a, i—c i. a—c 

= — (c\h\c) — V^(ci|ci) — (ci\ic) 
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=>- The eigenvalues of occupied Hartree-Fock orbitals are the negativ ionization energies, provided the 
orbitals are not allowed to change (relax). 

For the unoccupied states with orbital energy e r we can derive an analogous expression: 

EA = E(N)-E r (N+l) = -e r (3.55) 

Where EA is the electron affinity. 
Koopman's theorem: 



Given an N-electron Hartree-Fock single determinant ^fp with occuoied and unoccupied (virtual) 
spin orbital energies e a and e r , then the ioniszation potential to produce an (N-l)-electron single 
determinat W^-i w ith identical spin orbitals, obtained by removing an electron from spin orbital (p a 
and the electron affinity to produce an (N+l)-electron singe determinatn ^f/^x with identical spin 
orbitals obtained by adding an electron to spin orbital ip r are just — e a and — e r respectively. 



Table 3.6: Lowest ionization potentials (in eV) of some small 
molecules in HF. Values taken from A. Szabo and N. 
S. Ostlund, Modern Quantum Chemistry, Introduction 
to Advanced Electronic Structure Theory (Dover) 



ionization potential (eV) CH 4 NH 3 H 2 Q FH 



HF 


14.86 


11.65 


13.80 


17.69 


Exp. 


14.40 


10.88 


12.60 


15.81 


deviation 


3.1% 


7.1% 


9.5% 


11.9% 



Hartree Fock: 

• Overestimates ionization potentials (see tab. (3.2.5)) 

• Undersetimates electron affinities - often negative ions are not stable in Hartree-Fock 
In solids (see fig. (3.19): 

• Band gap severly overestimated, Hartree Fock ~ 6eV, experiment ~ 0.7eV 

• Shape of bands is correct, but bandwidth considerably overestimated 



3.3 Form of the exact wave function and configuration interaction 

In the previous section we had introduced a single slater determinant formed by occupying the lowest spin 
orbitals. But we can also occupy differently. 

f2k\ (2fc)! 

In total there are f ^ J = J^y^^ — ^pjj man y configurations. We have seen that a single slater determi- 
nant of Hartree-Fock spin orbitals gives the Hartree-Fock approximation to: 

How can we expand the full many-body wave-function ^e/ec? Suppose we have a complete basis {xi}- A 
function of a single variable can than be exactly expanded as: 

$(x 1 ) = ^a J x J (xi) (3.56) 



31 



3 Wave function based approaches 
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Figure 3.19: HF band structure of the Germanium crystal. Experimental values (squares) correspond to 
direct and inverse photoemission spectroscopy and are taken from Landolt-Bornstein, New 
Series, Group III, Vol. 17, Pt. a (Springer- Verlag, New York, 1982) and A. L. Wachs et al. 
Phys. Rev. B 32, 2326 (1985) 



A function of two variables can be expanded in stages: 

$(xi,x 2 ) = X a * (x 2 )xi (xi) 



(3.57) 



Inserting eq. (3.56): 



$(xi,x 2 ) = ^2hjXi (xi)xj( x 2) 



(3.58) 



If $ is supposed to be a fermionic wave-function we require antisymmetry: 

$(xi,x 2 ) = -$(x 2 ,xi) 



bij - bji ba — 



=>• $ (x 1: X 2 ) = ^ XI ^ ^ Xl ) X 3 ^ ~ Xt ^ X ? ( Xl )l 

i j>i 

= J2hjV2* SD (xi,x 2 ) 

An arbitrary antisymmetric function can be expanded in terms of all unique determinats formed in terms 

a complete set of one- variable funtions Xi ( x ) ■ 

This can be easily generalized to the N-body wafe-function: 

*e/ec = O),*0 + £ C :^+ £ C r X+ J] C :i t cKfc + --- (3-59) 

ra a<b, r<s a<b<c, r<s<t 
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X2k-1 



4 + 



Xik 



Xn+i -f- -|- Xn+2 

XN-l-\- -J-XiV 

X3 4~ -J- X4 

X2 4" 4" Xl 



Xn-i — -J- Xn+2 or Xr 
XiV+l-f- — Xn or Xa 

4 + 

singly exited 
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4- 



— |— Xs 

-J- Xr 

X6 

Xa 



4 ' + 

doubly excited 



*2 



Figure 3.20: Single slater determinats with different occupation. 
The coefficients are then obtained by diagonalizing. 



( V f elec\H elecl 1 ^ ' elec) 







S 


ab 

D 


aoc ' ' ' 

T 




(*o]fl"|*o) 


(*o}H\S) 


(*o\H\D) 




S 


(S\H\* ) 


(S\H\S) 


(S\H\D) 




D 


{D\H\9 ) 


(D\H\S) 


(D\H\D) 




T 


{T\H\* ) 









and the matrix elements are given by: 



(3.60) 







C*2 — X)i<j lri-r.fl 


case 1: 






\K) = \...mn...) 


(X|Oi|A0 = J2 m (m\h\m) 


(if 2 if) = \ Ysmn {mn\mn) - (mn\nm) 


case 2: 






\K) — | . . . ran . . .) 
|i) = |...pn...) 
(differ by one) 


(K\Oi\L) = (m|%) 


(K\0 2 \L) = J2 n (mn\pn) - (mn\np) 


case 3: 






\K) = \ ...mn...) 
\L) = \...pq...) 
(differ by two) 


(Jf|Oi|L) = 


(K\0 2 \L) = (mn\pq) - (mn\qp) 



and zero otherwise 

• The martrix elements are quite simple and involve at most the Coulomb integrals of 4 states. 

• Each slater determinant describes a configuration of electrons in the spin orbitals. Through the 
Hamiltonian these configurations interact with each other. 

=> Configuration Interaction (CI) 

For a complete basis CI would be exact but: 

• Basis is never complete 
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• (^) many configurations 

=>■ Complexity grows exponentially 
=> Only applicable to small systems 

• Defines correlation energy 

Ecorr = E e i ec — ^Eo^ (3.61) 

— Ehf 



HF limit exact result 




# of slater determinants 



Figure 3.21: Evolution of Hartree Fock and Correlation Interaction with respect to the number of basis 
functions and included slater determinants. The number of slater determinats can be reduced, 
according to the level of sophistication in the approach/approximation. 



3.4 A brief excursion into pertubation theory 

Another way to go beyond Hartree- Fock is to use pertubation theory. If we treat Hhf as zeroth order 
approximation to H e i ec we can write: 

H elec = H HF + V where V := H elec - H HF (3.62) 

or more generally: 

H e i ec =H Q + V (3.63) 
where Hq is any Hamiltonian we can solve easily: 

H ^=E^f or i/o|*! 0) )=^ (0) |*! 0) ) (3.64) 

For Hartree- FocK are the ground state slater determinants and all excited slater determinants formed 
from the Hartree-Fock single particle orbitals. If now is close to ^ e iec the effect of the pertubation V 
is small. 

=>■ Maybe there is an expansion that systematically improves the wave-function and the energy. 
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For this we introduce an expansion parameter A: 

H = H + XV 
and then expand in powers of A: 



A = H = H 

A = 1 H = H e i ec 



E. 



(o) 



\E, 



(i) 



A 2 K (2) + ... 



|^«=) = |*0) + A|*J) + A 2 |*?) + ... 



(3.65) 



(3.66) 
(3.67) 



We now wish to express the n-th order quantities in terms of the zeroth order quantities. For that we 
introduce the intermediate normalization: 



(3.68) 



= 1 and =0 V n > 

and insert eq. (3.66) and eq. (3.67) into eq. (3.65): 

{H + XV) (1*0) + A|tt}) + A 2 |v^2) + ...) = ^(0) + XE W + X 2 E m + ) + x ^i } + X 2^2 } + _) 

(3.69) 

Regrouping in terms of orders of A gives: 

n = /r |*?>=25j 0) |*?> 

n = 1 H \*l) + K|*°> = £7<°>|*J> + if } |*°> 

n = 2 H \^) + V\*l) = if } |* 2 ) + if^ 1 ) + if } |*°) 

Acting with (\l>?| from the left yields: 

if } =<*?|ff |*?) 
Ef> =<*?|V|*?) 

■^■Simple, because terms only involve \I>° 

eP =<*?|V|*J> 



For the ground state = *o and H = Hhf we get: 



HF 



E^ = (*? F \H HF \*n=j24 

i 

^->-mean field Hamiltonian acting on Slater determinant gives sum over eigenvalues 



4» = «>i*r> 



2 '—' |r, - r 



= i ^(mn|mn) — (mn\nm) — '^^(mn\mn) — (mn\nm) 



from cq. 3.20 

= — — ^^(mn\mn) — (mn\nm) 
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=> E HF = E^ + (3.70) 
For 2nd order we need an expression for |*o)- This we get from the expression for n = 1: 

(e^-h )\^\) = (v-e^)\^) 

= (V-<*?|V|*?))|*?> (3.71) 

This is no longer an eigenvalue problem, because no eigenvalue shows up on the right. We solve it by 
expanding |\I>|) in a basis. For this we take the solutions of H : 

n 

are obtained by acting with | from the left: 
From intermediate normalization we know that c-^ = 0. Reinserting into the expansion: 

i*?> = £*S><*M> ( 3 - 72 ) 

Multiplying equation (3.71) by gives: 

(e^ - E^) (^) = (^ n \V\^) (3.73) 
No we can insert eq. (3.72) and (3.73) into the expression for the 2nd order energy. 

(^\v\K)«\v\^) 



E 



p(o) F (o) 

E|(* |l/|*°)| 2 
l v llil 1 " ; l (3 74) 

P (0) p (o) ^-'^ 

The 2nd order energy is expressed in eigenfuntions >3/° and eigenvalues E^ of H , i represents the order 
of the excitation in the determinant. In other words: 

• Take the eigenfuntions and eigenvalues of H (the Slater determinant) 

• Form the matrix elements with the pertubations operator 

• Sum over expression (3.74) 

(2) 

If Hq is the Hartree-Fock Hamiltonian Hhf, Eq becomes (without proof): 



E {2) = 1 (ab\rs)(rs\ab) _ 1 (ab\rs) (rs\ba) 

2 ^ e a + e b -e r -e s 2 ^ e a + e b - e r - e s 

abrs abrs 



ab ab order ab ba order 

This is also called the M0ller-Plesset pertubation theory (MP2 energy). 
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In the previous chapter we have considered theories that work with the many-body wave function to solve 
the electronic Schrodingcr equation. However, the many-body wave function is quite an unwieldy object, 
with 3N-many coordinates. 

* = *({xi}) (4.1) 
Wouldn't it be nicer to work with objects that just depend on one variable like the density. 

n = n (r) 

This is not unreasonable because knowing the wave function implies we also know the density: 

n(r) = N J da J dx 2 ...dx N ^({xj})^ (4.2) 

normalized to 1 

J dm (r) = TV (gives particle number) (4-3) 

4.1 Kohn-Sham's equations 

Kohn and Sham's solution: 

=>■ Consider an auxiliary system of non-interacting electrons that has the same ground state density as 
the fully interacting system. 




Figure 4.1: Auxiliary system with same ground state as interacting system. 



h aux = -\+V{v,a) (4.4) 

The non- interacting electrons move in an effective potential V e ff that we assune to be V-representable, 
but this is a problem we know, because our full Hamiltonian 

H ({x 4 }) = ^2 h "-ux (xj) = H aux ({x 4 }) (4.5) 
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now becomes a sum over single particle Hamiltonians. These we know how to deal with from Hartree-Fock. 
hauxfi — tifi we get a set of single-particle orbitals 

=^dcnsity: 

N/2 

n ( r ) = E / dcr i^ ( r ' CT )i 2 

=^This defines the classic Coulomb interaction energy: 

E„H = \J^>"-^ (4.6) 
=^Thc kinetic energy is now also trivial: 

N/2 N/2 

Ts[n] = --J2 /d^|VV») = 2E l^fif ( 4 - 7 ) 

and 

£/rs N = T s [n] + J drV ext (r) n (r) +E H [n] + y«| (4.8) 

v Excf^l for now 

We now apply the same procedure as in the derivation of the Hartree and Hartree-Fock equations: 

L[n] := E KS [n] - E e t J dxtp* (x) tp t (x) 



,5L 



Vi 



+ 



SE^ + SEh + SExc 
Sn Sn Sn 



Sn 

— - enpi = 



^ - 2 V2< ^ 4 W [ v ext (r) + (r) + w X c (r)] (x) = (x) 



4 v2 + W0 



(r) = ei(fii (r) (only wave function depends on spin) 



The Kohn-Sham or effective potential is given by: 

Vks (r) = V eff (r) = v ext (r) + v H (r) + v X c (r) 



"ic (•") := 



5n (r) 



(4.9) 
(4.10) 



4.2 Hohenberg-Kohn Theorem 



Can we now recast 



Helec^elec = E e i ec ^ e i ec 
N 2 AT 



N v2 N 1^1 

E^r+E^* ^) + oEttz 



2 •'— ' |r, : - r,- 



ext 



+ Vet 



in terms of the density? It is in principle self-evident that the external potential, i.e. the positions of the 
nuclei, determines the properties of the system, but is this also true for the density? 
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Hohenberg-Kohn Theorem I 

The ground-state density n (r) uniquely determines the potential up to an arbitrary constant. 



Proof 

• The proof proceeds by contradiction 

• Let's assume non-degenerate ground states 

Suppose we have two external potentials v\ xt (r) and v 2 xt (r) that differ by more than a constant, but lead 
to the same density n (r) . 

v ext ( r ) T^ext ( r ) + const - 

H 1 H 2 two different Hamiltonians 

I 1 

Vl* 1 ty 2 two different wave- functions 

Then by the variational principle we have: 

(tf 2 |#l|*2) > = E 1 

as a side note: 

(*|T4 xt |*) = /dx 1 ...dx N ]Tiw(r)|*({x 1 ...x A r})| 2 

= ^2 drd<jdx 2 ■ • -dx N v ext (r) |* (r,CT,x 2 . . .xat)| 2 

i ^ 

= J drv ext (r) n (r) 

(* 2 |tfi|*2) = (*2|^ 2 |*2> + (* 2 |Hi - H 2 \^ 2 ) =E 2 +Jdr [v\ xt (r) - v 2 ext (r)] n (r) 
(T and V^ e give the same constant because the wave-funtion is the same) 
=> < E 2 + J dr [vl xt (r) - v 2 ext (r)] n (r) 
but we can swap indicies 1 and 2: 

=^E 2 <E 1+ J dr [v 2 ext (r) - (r)] n (r) 

=>■ There cannot be two different v ext that differ by more than a constant that give rise to the same 
density. 

=>■ This uniquely determines the external potential. 

Vext (r) <S=^ n (r) (4.11) 



Corrolary: 

Since the Hamiltonian is fully determined (except for a constant) the many-body wave-functions (for 
ground and excited states) are fully determined 

=>■ All properties of the system are fully determined given only the ground state density. 



39 



4 Density functional theory 



Hohenberg-Kohn Theorem II 

A universal functional for the energy E[n] in terms of the density n (r) can be defined, valid for any 
external potential v ex t (r). For any given v ex t (r), the exact ground state energy of the system is the 
global minimum of this funtional, and the density that minimizes the functional is the exact ground 
state density. 



Proof: 



1. Variational space =>■ restrict to u-representable densities (densities that can be represented by a 
potential v (r) 

2. Definition of functional => Density determines T and V ee 

EHK[n]= T[n] + V ee [n] + J drv ext (r)n(r) 
= F[n] + drv ext (r) n (r) 



Suppose m (r) is the ground state density of v\ xt (r) . 

=^.E l =E HK [ni] = 

Let's now consider a different density rii (r) that corresponds to a different wave-function 4 , 2- 

E x = E HK [ ni ] = < (* 2 |ffi|*a> = ^2 

Hohenberg-Kohn functional evaluated at ground state density gives the lowest energy 

=> If the functional of the density is known, than the total energy of the system can be obtained by 
variational minimization with respect to the density 

==> Hohenberg-Kohn functional only gives ground state, not excited states (like e.g. Configuration Inter- 
action) 



Hohenberg-Kohn theorem 
establishes one-to-one correspondence 
between v ext and n 

Vext (r) 



Determines H 
and gives all 
many-body wfn. ^ f 



We pick the 
ground state 



n (r) 
A 



That determines the 
ground state density 



Figure 4.2: Short Summary 

Problems: 

• Our proofs somewhat want back to the many-body wave-function 

• There exists no prescription to determine the kinetic energy from the density 

• We have no prescription for generating densities 
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EhkH = F[n] + J drv ext (r) n (r) (4.12) 

Eks M = T s [n] + J drv ext (r) n (r) + E H (r) + E xc (r) (4.13) 
^E xc \n] = T[n] - T s [n] + V ee [n] - E H [n] 

1. 2. 

1. Difference between interacting and non-interacting kinetic energy 

2. Difference between electron-electron interaction and classic Coulomb energy 

• The key now is to find good approximations for ExcV 1 } 

• Exc [n] is in general small, because T$ and Eh capture a large part of T an V ee 

Simple approximations might already be successful 
! Altough what matters again are energy differences and there Exc [n] can be decisive 

4.3 The local density approximation 



n (r) . 




Figure 4.3: Slowly varying density. 

For inhomogeneous systems that have a slowly varying density the system localy looks like the homo- 
geneous electron gas (see fig. (4.3)). 
==> constant density n constant external potential 



E xc [n] = / dm(r)ea; C ([»],r) 



(4.14) 



e xc is the energy density (i.e. the energy per electron) at point r in space, that depends only on the density 
at that point. The exchange-correlation potential v xc follows straight forwardly by functional derivative: 



, lc (r).f#. £ (Hr) + n(r) ie([lil ' r) 



Sn (r) 



Sn (r) 



(4.15) 



The exchange-correlation enrgy density for the homogeneous electron gas is known from accurate Quatum 
Monte Carlo calculations: 



e x n 



3kp 

47T 



1 /3 

with the Fermi momentum hp = (37r 2 n) 



(4.16) 
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is given exactly. 



and for v T 



E xc [n] = J drn(r)e xc ([n],r) = / drn(r) e IC (n(r)) 
= J dm(r)|-(3^ 2 n) 1/3 



3 /3 



4 V 7T 



3 /3 



Very simple expression! 



1/3 



1/3 



dm (r) 



4/3 



,1/3 



3 / 3 \ 1/3 1^-2/3 



43 / 3 
34 ^ 

3 x 1/3 
7T 



1/3 



,1/3 



,1/3 



4 6 



r = (_3_\ 1/3 = 1-919 

s \47rn/ kf 



Figure 4.4: The correlation energy density of the homogenous electron gas. 

The correlation energy density for the homogeneous electron gas is not known analytically, but can be 
computed to very high precision using Quantum Monte Carlo techniques (like Configuration Intcration 
it is a method that works directly with the many-body wave- function). An analytic expression for this 
behaviour (see fig. (4.4)) was first estimated by Wigner in 1938: 



0.44 



r s + 7, 

later better parameterizations by e.g. Perdew and Zunger: 



(4.17) 



Gcll-Mann/Briickncr 



^41nr, 



B + Cr s \nr s 

7 



1 + piJrS + fcr. 



r s < 1 
r, > 1 



(4.18) 
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(The numerical values of A, B, C, D, 7, fa and fa are given in Perdew, Zunger, Phys. Rev. B 23, 5048 
(1981)) 

• Spin polarized extension 



E 



or if we introduce the spin polarization: 



^ [„t „4-] = |drn(r)e xc (nW) 



(4.19) 



£(r) 



rv (r) — (r) 
ra(r) 



£f c iz?A [nt, n+] = / drn (r) e xc (n (r) , £ (r)) 



(4.20) 

What we need now is reference data for the spin-polarized homogenous electron gas and appropriate 

paramterizations, which can also be found in Perdew, Zunger, Phys. Rev. B 23, 5048 (1981). 

LDA: 

• Is by construction exact for HEG 

==> Expected to perform well for systems with slowly varying density (e.g. simple metals) 

• Typically gives dissociation energies of molecules and cohesive energies of solids within 10-20% 

• Bond length (lattice constant) typically within 1-2% (and typically to small) 

• Problems e.g. for rapidly varying densities e.g. atoms 



4.4 Generalized Gradient Approximation 



Ar Atom 



5~ 




Figure 4.5: Radial density in Ar Atom, From The ABC of Density Functional Theory by Kieron Burke 



• The shell structure is clearly visible in figure (4.5) 

• The density is far from homogeneous! 
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The idea: add gradients to E xc : 

E xc A [n 1 ",^] = J drn(r)e xc (n^n 4 -, |Vn f | , |V^| , . . .) 

= J dm (r) e% EG F xc (V,^, |Vn T | , |Vr^| , . 



Generalized gradient expansion: 

It makes sense to work with a scaled gradient 



S(t) 



|Vn (r) 



2fci?n (r) 

that measures the gradient on the scale of the density itself. 



Eg GA [n\ n*>] = J dm (r) e H EG [n]F xc (n\ n\ S\ 



(4.21) 



(4.22) 



(4.23) 



03 




Figure 4.6: Scaled gradient for Ar Atom, From The ABC of Density Functional Theory by Kieron Burke 

• There's no unique form or parameterization for F xc 

=>■ Many different parameterizations have been proposed by now that fall largely into two cate- 
gories: 

— Satisfy a certain number of exact constraints 

— Are fit to benchmark sets 

• GGAs usually work best in the regime for which they were designed 



4.5 Self-interaction and exact-exchange 

Our DFT energy expression is: 

E to t[n] = T s [n] + E ext [n] + E H [n] + E xc [n] 
Let's recall that in Eh we summed over all single particle states in the system for convenience: 

_,#(r)V< (r)^(rOVi(r') 



(4.24) 



^ r i 1 f , , t n (r) n (r') 1 
E H n] = - / drdr' y V = - 
1 J 2./ r-r' 2 



? / dl " dr ' 



r - r' 



= i£/drdr< 



,l^(r)| 2 |^(r')| 2 
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so when i = j an electron interacts with itself =>■ self-interaction 
In Hartree-Fock this term is exacly canceld by the exchange energy: 



E.. 



drdr 



>*(r)^(r)V*(r')^ (r') 



For i 7^ j we now have terms 



2 J |r - r' 

that exactly cancel those coming from However, for LDAs and GGAs this cancellation does not occur: 



E LDA/GGA = J drn (r) ^ ( ^ r (Vrj _ _ _ )} 

= y dr|Vi (r)| 2 £;cc ([n],r,(Vr,...)) 



Perdew and Zunger (in their 1981 LDA paper) defined a 1-electron self-interaction error on this basis: 



2/ r-r' 



Ilk (r) 



(4.25) 



• <5j is in general not zero 

• As we saw earlier in section (3.1) the self-interaction error has a tendency to delocalize states 

We know that Hartree-Fock removes the self-interaction error but can it be cast into the DFT framework? 
For that we would need the functional derivative with respect to the density: 

? HF 



E = T s + E ext + E H + E x 



The challenge is 



all other terms arc standard, but 



chain rule: 



and chain rule applied again: 



v x (r) = 



SE X 



Sn (r) ' 
n(r)= J>< (r)| 2 

i 

v fr) = Wdr' SE ' 

v x (t) # . (r0 Sn(r) 



+ c.c. 



(r)=J2 J dr ' dr " 



SE X 5jH (rQ 
Stpi (r') 8v KS (r") 



8v KS (r") 
Sn (r) 



The first term we have already encountered in the derivation of Hartree-Fock: 
6E X ™ f , >*(r')^(r')^(r') 



<lr' - - - ^^^ - / dr% (r, r') ^ (r') 



The second term 



Hi (r') 

Sv KS ( r //) 



(4.26) 



(4.27) 



(4.28) 
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follows from first order pertubation theory. We recall: 

e l = e° + \e^ + ... 
4>i = + A^ (1) + • • • 
vks = v° + Aw (1) + . . . 



inserting this in the Kohn-Sham equtions 



h KS ipi = aipi with h 



KS 



vks 



gives: 



n = h° KS ^ = e°V° 

n=l h° KS ^+v^ = e^+e^^ 



now we need to get rid of e 



(i) 



>apply (ip°\ 

^^l^fll^ (1) )+<^l« (1) l^ 0) )=e?<^l^ (1) > 



Now we recall the expansion of ip^ in terms of {ipj}'- 



m = e - E ( yf if) 



From this follows with ch/>i = IV^) an d <5v/fs = 

tyi (r)_ ^ V.* (r') ^ (r) 

The last term 



r ') & ^ 



^ (r') = Gi (r,r')^ (r') 



5v (r) 



<Jn (r') 



is the inverse of the Kohn-Sham response function: 

|2 N 

Sv (r') ^-r J 5v (r') 

JV 



^ = E^r=E^(rO^(r, r ^ l ( r O + e.c. 



EE + c - c - 



= X S (r,r') 

Putting this all together yields: 

vx (r) - E / dr ' / dr " / dr "' E * ( r '' r ") ^ ( r "') °i ( r '< r ") ^ ( r ') + c ' c - 
This expression 



(4.29) 



(4.30) 
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• Is rather involved and computaionally challenging 

• But it shows that a self-interaction free multiplicative potential can be constructed 

This formalism is also known as OEP: optimized effective potentail approach, because v x is the variationally 
best potential to E x or E xx in solids 

• E x on\y implicitly depends on the density 

E x =E x [{^ i }]=E x [{^ i [n}}] 

This is known as orbital functional 
OEPx or EXX (see figure (4.7) and (4.8): 

• Moves Hartree-Fock into the realm of DFT 

• Is self-interaction free 

• But includes no correlation 



— — 4p ... — 1 

3d 4s ... 5p z 

57 j 

■■■ LDA 2p 

EXX -_ 

2s~ j 

Scandium Nitrogen Indium ; 

Figure 4.7: Kohn-Sham eigenvalues of the highest atomic states in LDA and EXX (=OEPx). (from Qtcish 
et al. Phys. Rev. B 74, 245208 (2006)) 



> 
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4.6 Hybrid functionals 

OEPx (EXX) and HF contains too much exact echange. A pragmatic solution that was adopted first in 
quantum chemistry was to construct a hybrid: 

E h x f = E° C FT + a (E» F - E? FT ) (4.31) 

In this simplest form a portion of DFT exchange is replaced by exact exchange, while correlations remain 
on the DFT level. More complex hybrids with more parameters exist (e.g. B3LYP) or where E^ F is range 
restricted (long, short). One of the most popular choices is a — 0.25. Combined with PBE-exchange this 
functional is known as PBE0 (Perdcw-Burkc-Enzcrhof). Contrary to common believe hybrid funtionals 
are not cast into the Kohn-Sham formalism of multiplicative potentials via the OEP formalism. Instead 
^fffi 2 is performed like in HF theory, which leads to a non-local potential: 

(r, r') = [v^ T (r) - av D x FT (r)] S (r - r') + a^ x (r, r') (4.32) 
The non-local potential has certain advantages over the multiplicative potential, as we will see later. 
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Figure 4.8: Kohn-Sham band structure of ScN and InN. In LDA ScN is a semimetal and InN a metal. 

EXX correctly predicts both compounds to be semiconducting, (from Rinke et al. phys. stat. 
sol. (b) 245, 929 (2008)) 



4.7 DFT and excitation energies 

In general, and for excitation energies in particular, it is very important to distinguish two things: 

• Can a quantity in principle be calculated exactly with exact DFT? 

• How do approximate funtionals perform? 

Let us for now consider two different simple excitations (figure (4.9)): 



cxitation energy 



e s = E 



tot 
final 



E 



minimal 



(4.33) 



1. e s = E(N,s) - E(N) 

• For a wave-function based methos (CI, Coupled Clousters, etc.) es are the eigenvalue differences 
of = 

2. £5 = E(N — 1, s) — E(N) - ionization energies 

3. es = E(N, s) — E(N + 1, s) - affinities (are for some reason defined as Ei — Ef) 
2. und 3. do not come out of = E^. 

Where are we with our two questions? The way the excitiation energies are written they involve only 
differences in total energies. But remember in DFT only the ground state can be exact. 

only 

• ionization potential I = E(N — I) — E(N) 

• electron affinity A = E(N) - E(N + 1) 

• gap E gap = I - A 
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Figure 4.9: Simple excitations 



could be calculated exactly in exact DFT. 

All quantities involving an "s" are excited states of the system, that cannot be expressed as difference 
of two ground states. To mention only two problems with trying to prepare excited states E(N, s) or 
E(N±l,s): 

• Finding a suitable constrained to keep the system in state s may not be possible or it may not survive 
the self-consistency cycle 

• Excited state densities are not unique (unlike ground state densities), i.e. there is no Hohenberg- 
Kohn theorem for excited states 

What about the Kohn-Sham eigenvalues then? In Hartree-Fock we had Koopmans' theorem: 

-e* F = E*{N±l,s)- E{N) 

Where the orbitals for the calcultaion of E* were constrained to be the ones of E(N) (frozen orbitals). In 
DFT there is no such theorem/relation. It can, however, be proven that the highest KS eigenvalue of a 
finite system equals the negative of the ionization potential. 

1 = -£n S ( n ) in exact KS (4.34) 

A hand-waving argument for this is: 

• The asymptotic long-range density of a bound system is governed by the occupied state with highest 
eigenvalue 

• Since the density is supposed to be exact, so must the eigenvalue be 

(A more rigorous proof can be found in C. Almbladh and U. von Barth, Phys. Rev. B 31, 3231 (1985) or 
M. Levy, J. P. Perdew and V. Sham, Phys. Rev. A 30, 2745 (1984)) 



For approximate functionals this is not true, see figure (4.10) 
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Ionisation potential in the LDA 



0.0 
-10.0 
E£ -20.0 
£ -30.0 
-40.0 
-50.0 
-60.0 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 










#-# E(N-1)-E(N) 






O-O -e N ( N ) 













He Be Mg Ca Sr Cu Ag Au Li Na K Rb Cs Zn 



Figure 4.10: Ionization potentials of atoms calculated as total energy difference in LDA and by the LDA 
Kohn-Sham eigenvalue. Shown is the error with respect to the experimental ionization po- 
tential. (Reference: NIST - Atomic reference data) 



Then there is Janak's theorem that establishes a connection between the KS-eigenvalues and the deriva- 
tive of the total energy: 

B E 

— = e, Janak, Phs. Rev. A 18, 7165 (1978)) (4.35) 
where n is the occupation of a given state i. For the proof we write: 

U = J drf* (~\^ = 6i ~ J dr ^i ( Vext + V H+ Vxc) *i 



JV 



drii ^— ' drii J ^— ' drii 

3 3 » ^~ 



fi, ' 



Now considering the first part of eq. (4.36): 



dtj f 0*J / V\. 



(4.36) 



Then we introduce the occupation factors nc 

N N 

n{r) =^n i |* i | 2 and T^^mU (4.37) 

i i 

^E = f + E ext [n] + E H [n] + E xc [n] (4.38) 
^L = ti +J2nj^: + J dr(v ext +v H +v xc ) |^| 2 +^ 

3 \ 3 

Inserting the second part of eq. (4.36) into this equation cancels the first part in the integral. 

()E x - dtj f x - ' ''' 
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dE r 



5** /V 

"7^- I + Vext +V H + V xc ) Wj +C.C 



5** 

— + c.c 



9 f 

£, + nj£j — — / dr | SI/ j | 2 • 2 (wave- functions are normalized) 
on>i J 



=o 



= c< □ 

Rearanging Janak's theorem: 



E(N + — E(N) = / dnei(n) ps e(0.5) (mid-point approximation) 
j o 



In other words: 



(4.39) 



• Excitation energies are approximately given by the value of the eigenvalue of half occuption (also 
known as Slater- Janak's or Slater transition state) 

• The problem remains, however, that for any but the highest (lowest) state occupations need to be 
suitably constrained. 

Derivative discontinuity 

Let's consider the gap of a large system (finite, but so large it could be a solid). The gap is defined as: 




N electrons 



N+ 1 electrons k 



Figure 4.11: After the addition of an electron into the conduction band (right) the xc potential and the 
whole band-structure shift up by a quantity A xc . Figure after R. W. Godby et a/., in A 
Primer in DFT, Springer 2003 



E gap =I-A = E(N + 1) - 2E(N) + E(N - 1) 
But we know that the highest occupied state in exact KS is exact: 

rp J<S { at i 1\ ^KS 

hi. 



(4.40) 



= e« s +1 (N + 1) - e§ s +1 (N) + e^Af) - e 



E KS 
gap 
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Because our system is large we have 

N > 1 => An(r)^0 for N -> A + 1 

In other words, our density change is infitesimal. 

=>■ ujj and u e2:t will not change 

=> A xc can only come from change in v xc 

=>■ the derivative of E xc with respect to particle number changes discontinously 
=> v xc changes by a constant (see figure (4.11, 4.12)) 




r(a.u.) 

Figure 4.12: Exact-exchange potential (KLI) shifts up almost uniformily by w 1 eV upon adding 10 -6 
electrons. J. B. Krieger et al, Phys. Rev. A 45, 101 (1992) 

• Even exact KS calculation will not capture the derivative discontinuity and therefore will give incor- 
rect gaps, if only KS eigenvalue differences are considered 

• There is currently fierce debate about the size of the derivative discontinuity in exact KS 

• The role of the kinetic energy is unclear 
In summary: 

• We have certain exact relations linking total energies and eigenvalues with excitation energies 

• But then there is the derivative discontinuity and approximate functionals (LDAs and GGAs do not 
contain the discontinuity, a problem even for total energy differences and they are plagued by the 
self-interaction error) 

• Note that non-local functionals do not really have the discontinuity problem because they are not 
generated from |f (see figure (4.7, 4.8, 4.10, 4.13)) 
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Figure 4.13: Typical behaviour of most current DFT functionals for the case of a finite system (e.g. a 
cluster) that approaches its bulk limit. 
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Figure 5.1: Time dependence of excitation process 

• Include time dependence of excitation process 

• Switch to time dependent pertubation theory 

• What about total energies? 

5.1 Second Quantization 

In this section we introduce the method of second-quantisation, which provides an elegant, alternative to 
the conventional notation in quantum mechanics. In second-quantisation a physical state is represented by 
a state vector in Hilbert space and observables by an operator acting on this state. If the operator applied 
to a state reproduces the state bar a multiplicative factor, we speak of an eigenstate. The eigenstates of 
a Hermitian operator form a complete set. With this knowledge we can chose a suitable, complete set of 
states to represent the many-body problem. In the set of all Slater determinants based on a complete set 
of single particle functions ipk(x) we find a basis system and a many-body wavefunction that meet the 
requirements. We label this many-body state and use the short hand notation 

$N = \ki,k 2 ,...,k N > (5.1) 

The fermionic creation operator a\ creates a particle in state k if this has previously been unoccupied 

a\\ki,k2, ■ ■ • ,k N >= \k,k-L,k 2 ,...,kN > (5.2) 
similarly a particle is destroyed in state k by the annihilation operators a k 

a k \k,k ll k 2 ,...,k N >= \ki, k 2 , ■ ■ ■ , fcjv > (5.3) 
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Based on the definitions for a\ and a k we define the field operators tp^{x) and V(x) 1 

^t(x)=^4^(x) (5.4) 

k 

^( X ) = ^2a k tp k (x) (5.5) 

k 

The anti commutation relations for the creation (5.2) and annihilation (5.3) operators are retained for the 
field operators: 

{*(x),^(x')} = S(x,x') A {V t (x),V t (x')} = {V(x),V(x / )} = (5.6) 
In field operator notation operators of the form Oi(x) and C*2(x') assume the shape 

Oi=J2 f rfri/- t (x)0(x)^(x) (5.7) 

cr •* 

2 = J2 I ^r , 7/; t (x)V t (x')0(x,x')V(x)V'(x') (5.8) 
The Hamiltonian in second quantized form 

H = Y^(i\h\j)a\aj + \ |W)at a ta fc ai (5.9) 

ij ijkl 

is thus transformed according to 

H = E/ dr^( X )h(r)i>( X ) + \J2J J dvdv' ^\x)^(x!)y^-^{x)^) (5.10) 

where h(r) is the single-particle Hamiltonian 

h(r) = -^V 2 + v ext (r) (5.11) 

as always. 

For the following analysis it is convenient to switch to the Heisenberg notation for the field operators 

V>(x,i) = e tHt %b(x)e- tHt (5.12) 

If we now consider the iV-electron ground-state, or in short \N > and an eigenstate s of the N + 1 
particle state: | N + 1, s >, we then obtain from (5.12) 

(N\1>(x,t)\N + l,s) = (N\e lHt iP(x)e~ tHt \N + l,s) = e^'j^x) (5.13) 

with the excitation energy e s and amplitudes / s (x) defined by 

e s =E(N + l,s)-E(N), /s(x) = (N\^(y)\N + 1, s) (5.14) 

From the definition of the field operators (5.4,5.5) combined with the relation (5.3) we find that / s (x) 

W(x)|iV + l,s) = (N\ OfcVfc(x)|iV + 1, s) = (N\ V s(x)\N)8 sk = p a (x) (5.15) 

k 

gives the excited state wavefunction. In the next section we explore the connection between </? s (x) and 
the propagation of a quasiparticle is made. 



1 We introduce the shorthand notation x to denote the set (r, cr) 
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G 



t > t> © 



x',t' 




Figure 5.2: Single Particle Greensfvmction for t > t' and t < t' . 



5.2 Green's function theory 

5.2.1 The Single-Particle Green's Function 

We will begin our introduction by denning the single particle Green's function, G, and then proceed to 
show that the one-particle excitation spectra is naturally contained in G. 

Using the notation of second quantisation we define the single-particle Green's function as: 

G(xi,xY) = -i{N\f{ip{x,t)^ ,1/)}\N) (5.16) 

where ifv (x' , t') and '(/'(X; t) are the creation and annihilation operators of a particle, respectively, and the 
time ordering operator, T, assures an ascending order in the set of field operators from right to left. Every 
pair commutation of fermionic field in Equation (5.16) operators is accompanied by a sign change and 
(5.16) can thus also we written as 

G(xt,xY) = - i(N\ip(-x,t)^ (*.' ,t')\N)G(t - t') 

+ i(N\ip* (x' ,t')^(x,t)\N)e(t / - t) (5.17) 

For times t > t' the Green's function describes the creation of an additional electron at x' and subsequent 
propagation and annihilation at x and time t. Conversely for t < if a hole is created by extraction of an 
electron from the ground-state by the operation ip(x,t)\N >. The overlap with ip'(x.')\N > then gives the 
probability of the hole having propagated to x' at time t' . 

Making use of the Heisenberg notation (5.12) equation (5.17) becomes 

G(xt,xY) = - i(N\e lHt iP(x.)e- lHt e lHt ' ^ (x.')e~ tHt ' \N)0(t - t') 
+ i(N\e lHt V (x')e- lHt ' 'e tHt ^)e~ lHt \N)e(t' - t) 

= - ^(A^V>(x)e~ ^(H ~ £N)(^ ~ t 'V t (x')I^O ( ^(^ - 1') 

+ i(jV|V>V)e i(ff ~ B " )(t_t 'V(x)|i\Oe(i' - t) (5.18) 

If we insert the completeness relation in Fock-space 

1 = \vac X vac\ + £ \* l s X *J| + . . . £ |*f X *f | + . . . (5.19) 

s s 

where \^>f > denote the eigenfunctions of a N-particle system \N,s >, equation (5.18) transforms to 
G(xi,xY) =- iJ2( N ~ l\Tp(*)\N,s)e-^ EN - 1 - EN K t ~ t ">(N,s\iPHx.')\N ~l)e(t-t') 

S 

+ iJ2( N + 1|^V)I^ s )e l ( E »+i- E ^ t ~ t ">(N, s\1>(x)\N + l)6(t' - t) (5.20) 

s 

The sum over particle number disappears because the scalar product of wavefunctions with different 
particle numbers vanishes. In a more compact notation the Green's function appears as 

G(xt, xY) = -i A(x)/ S * (x')e- i£ « [0(t - t')Q(e s — /x) — Q(t' - t)Q(n - e s )] (5.21) 
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where we have defined the excitation energies e s and wave functions f s (x) by 2 



e s = E{N+l,s)-E(N), /s(x) = (N\i/>{x)\N + 1, s) for e s > M 
e s = E(N) — E(N — 1, s), /s(x) = (AT — 1, s\ip(x)\N) for e s < M 



(5.22) 
(5.23) 



The energies correspond to the true single-particle excitation energies upon addition and removal of an 
electron from the system in the ground-state. 

5.2.2 Spectral Representation of the Green's Function 

In order to extract more physical information from the Green's function we switch to the spectral, or 
Lehmann, representation of G. We first note that for explicitly time-independent Hamiltonians the Green's 
function only dependents on the time difference r — t — t' . Secondly wc need to treat the discontinuity in 
the time argument, brought in by the Heavyside function, carefully We state without proof that 6(±t) 
can be represented by 



6(±r) = lim =F- — 
v ' 2m 



duj 



uj ±in 



-i(e s +ui)r 



9(e s - n) - e s ) 



Inserting this equation into (5.21) we obtain 

/OO 7 
-oo 2m [ uj + if] uj - ir] 

and with the additional substitution uj = e s + uj wc identify the Fourier transform of G to be 

9(e s - (j,) e(fi - e. 



G(x,x»= lim T/ S (x)/J(x') 

7J— >0+ "—^ 



uj - (e s - ir\) uj - (e s + in) 



(5.24) 



(5.25) 



(5.26) 



II 



Poles of G 



XXXXXXXXXXXXXXXXXX) 



Im uj 



xxxxxxxxxxxx 



Re uj 



flN-l A*JV 



III 



IV 



Figure 5.3: The poles of the Green's function for a finite system lie infitesimally close above the real axis for 
energies smaller than the chemical potential of the N-l-particle system /ijv-i. For energies greater 
than /ijv the poles fall infitesimally close below the real frequency axis. In the limit of zero temperature 
a single chemical potential can be defined by /i, as indicated in this figure. 

From the spectral representation (5.26) wc deduce, that the poles of the Green's function are the exact 
excitation energies (5.22) and (5.23) of the system, referenced to the chemical potential. For a finite sys- 
tem the energy spectrum is discrete with the poles of G lying infitesimally above the real frequency axis 
for energies smaller than the chemical potential and infitesimally below otherwise, as illustrated in Figure 

2 In principle ft has to be defined for the addition (At(jv)=B(JV+l)-B(JV) ) an< J the removal (h(n-1)=e(N)-E(N-i)) of an 
electron separately. We will comment on this point in Section 5.2.2 



57 



5 Green 's function theory 



5.2.2. The Green's function is analytic in quadrant 3 I of the complex plane for energies larger than jitjv-i 
and in quadrant III for energies lower than /z/y For an infinite system A w A — 1 and the two chemical 
potentials merge into one. In this case the discrete series of poles form a branch cut in the complex plane, 
except possibly for a band gap region where there are no eigenstates. 



5.2.3 Expectation value of single particle operators 

According to equation (5.7) we have: 

Oi=J2 [ dri/. t (x)0(x)^(x) 
Let us write 0\ with non-locality in the spin index for now: 

oi = Y,J dr 4 ( r )°^' ( r )^« ( r ) 

af) 

the expectation value is then 

(iV|0 1 |iV) = (iV|^ f dr^(r)0^(r)^ a (r)|iV) (5.27) 

we now have to: 

• Pull the operator out of the expectation value 

• Swap the order of the field operators 

• Introduce the time ordering operator 
Using equation (5.4) and (5.5) we can write 

(N\tf (x) 6 (x) V (x) | TV) = (N\tf ( X ) 6 (x) J2 <Pi (x) aj\N) 

3 

We introduce an artificial x' dependence 

V> f (x) = lim ^ (x') (5.28) 

x'— >x 

==> now we can swap O and tp^ 

(N\tf (x) 6 (x) V (x) | A) = f dr Km ]T O af} (r) (A|^ (r) 4> a (r) | A) 

Then we can employ the limit trick once more: 

1 = lim e- 1 ^*-*') = lim c lH{t '-^ (5.29) 
t'->t t'->-t 



3 To be absolutely precise, it is only sensible to speak of quadrants in connection to the Green's function or the self-energy, 
if the chemical potential coincides with the origin of the complex plane. In the following Chapters the chemical potential 
is implicitly assumed to be zero unless otherwise stated. 
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(N\ipt (r') Va (r) \N) = lim e^*"*') (N\1>1 (r') lim e^*'"*^ (r) \N) 

= lim(N\e iH ^ipl {Y r )e- iH{t '^c lH{t ^ a {r) C -' lH ^\N) 
t'->t p 

= lim(7V|^ (r',t')i; a (r,t)\N) 



t'-n 

= lim(-l)(iV|Va (r,i) Vj, (r',t')l^> 



= Um(-1)<JV|T^ Q (r,i)^ (r'.OM 
t'>t 



- lim iG Q(3 (r,t; r',i') 



Putting this together: 



(AT|6|A0 = -i / dr lim lim V O Q0 (r) G Q ^ (r, i; r', i') (5.30) 

a/3 

-i y dr lim tr{o(r)G(r,t; r',t')} (5.31) 



The expectation value of any one particle operator can be determined with the one-particle Green's 
function. 



In particular: 

r v 2 

(T) =+i / dr lim — tr{G(r,t; r',t')} (5.32) 

J r'-s-r 2 

and for the density we get 

O a p(T)=8 aP 5(T-i f ) (5.33) 
n (r) = (7V|0|iV) = -i / dr lim V 5 afj 5 (r - r') G af} (r, i; r', i') 

/ r' — >r 

J a/3 

= -iG aa (r, i; r, t') = itr {G (r, t; r, t')} (5.34) 

what about the two particle operator, i.e. the Coulomb potential and therefore the total energy? 
For this we consider the equation of motion for the field operator: 

iJ^(x,t) = foHx,t),tf] (5.35) 

The commutator is calculated in the Heisenberg picture ((5.12) and (5.10)). Applying the anti commuta- 
tion relations (5.6) and the identity 

[A, BC] = {A, B}C — B {C, A} (5.36) 
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we obtain for 



(x),h} = J <^[y>(x),y>ty)fe(^y(x)] 



/ 



{^(x),V t (x')}/ i (x / )V(x / )-V' t (x / ){Mx / )V'(x / ),V'(x)} 



V 



(5(x-x') 



= /i(x)V(x) - /" dxV f (x') (/i(x')V>(x')V(x) + V(x)/i(x')V(x')) 

= /i(x)V(x)- /* dx' ((-l)V t (x')^(x / )V'(x / )V'(x) +^ t (x')^(x)/i(x')V(x')) 



= /i(x)V(x)- / dx' {V t (x'),V(x)}/i(x / )V(x / ) 



= /i(x)V>(x) (5.37) 
For the two particle part of id we obtain an analogous expression. The equation of motion therefore reads 



fc(r) + y V t (x',t)w(r,r')V'(x',t)dx' 



V>(x,t) 



(5.38) 



with u(r,r') being the Coulomb potential. 

We now multiply from the left with ?/>t(x',f') and take the ground state expectation value: 

(i9 t -/ i (x))(/V|V t (x / ,t / )V(x,t)|iV)= / dx"(N\^,t')^(x",t)V(v,v'))^,t))^,t)\N) (5.39) 



integrating over / dx gives the expectation value of the two particle 
part of the Hamiltonian in the limits x' — > x and t' — > t + 

(N\v\N)=-z I dx lim lim ( ^ - h (x) ] G (xt, xY) 

2 J t'— >*+ x'— 5-X \Ot J 

where we applied the same tricks as before to turn (N\ip^(x.', t')^>(x, t)\N) into the Green's function. Now 
we have: 

E a = (N\h + v\N) (5.40) 

recall that 



(N\h\N) = ( dx lim lim h(x)G {xt, xY) 

J t'->t+x'->x 



-i/ 



dx lim 

x'— >-x 



G , (xi,x't+) 



The ground state total energy can also be calculated from the single particle Green's function! 

Summary of this part: 

The one-particle Green's function gives: 

• Charged excitations 

• The ground state total energy 

What remains to be done is to find suitable approximations. 
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5.3 The self-energy 

Before we proceed to find approximations it is useful to introduce the concept of the self-energy. Going 
back to the equation of motion for the field operator (eq. (5.39)). 

(id t -h(x))(N\^(x',t'W(x,t)\N) = J dx"(iV^(x',^ t (x,t)F(r,r'))V(x',t))V(x,i)l^> 
We can apply similar tricks as before to write this as: 

(id t - h (x)) G (xi, xY) = 8 (t - t') 5 (x - x')-i / dx"w (r, r") (N\T (x', i)V(x", i)V(x, i))^ f (x, i')} |iV> 

J „ ' 

G 2 (xt,x"t,x"t+,x't') 

The two particle Green's function defined as: 

G(xi,x 2 t 2 ,x'i',x' 2 i' 2 ) = (AT|T{V(x 1 ,t 1 )V(x 2 ,t 2 )V t (x 2 ,t 2 ))V t (x / 1 ,t / 1 )} |iV) (5.41) 

Therefore, to calculate the one-particle Green's function requires knowledge of the two-particle green's 
function, which in turn will require knowledge of the three-particle Green's function and so on, It is easily 
seen that this builds up the full many-body Schrddingcr equation. Alternatively we can here make an 
attempt at factorizing: 

-i J dx%(r,r")G 2 (xi,x"i,x"i+,x'i') = J di"dx"M (xi, x"t") G (x"t", x'i') 

From M we seperate out the Hatree potential. 

v H (r)= J dr'v (r, r') n (r') = J dr'v {r,r') (N\tp^ (r' ,t')il>(r' ,t)\N) (5.42) 

So that with 

E = M + w H (5.43) 

we arrive at 

^ - h (r) - jiff (r)^ G (xi, x'i') = 5 (t - i') 5 (x - x') + J di"dx"E (xi, x"i") G (x"i", x'i') (5.44) 

and Fourier transformed: 

(w - /i (r) - w ff (r)) G (x, x', w) - J dx"E (x, x", w) G (x", x', w) = J (x - x') (5.45) 
If E (x, x",w) were static like in HF or local as in DFT this equation would simplify: 

(w-/io)G(x,x',w) = £(x-x') where ft - ft (r) + v H (r) + J VHF j*' ' 'J (5.46) 



"xc (r) <5 (r - r') 



AT 



H = 5>o (i) 



This is the homogenous differential equation to the Schrodinger equation 

ff *s - ^ 0) *. hofs = e a ip, (5.47) 
For this system we know that \N) is a single slater determinant 4>q. 

=^/ s (x) = (AT-l, s |V(x)|$ ) 

= (N- l,s\^2<f k (x)a fe |$ > 
k 

= VVfc (x) (N - 1, s|a fe |$ ) = *Ps (x) 

k 7 
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Go (x,x»=limy^ 



Lf s (x) ip* s (x') 



»7->o ^ w - (e s ± irj) 



for occupied 
for unoccupied 



states 



(5.48) 



So we have 



(w - /*,)G = J -»• (w - fto) = Gq 1 
(w - /i )G - EG = J 
rv G^G - EG = / 



G — Gq GqEG 



This is Dyson's equation! 

Starting from a non-interacting Green's function, G can be obtained through Dyson's equation. 

5.4 Hedin's equations 

Without proof Hedin's equations are listed here. The derivation can be found in the appendix of Hedin's 
original paper: Phys. Rev. 139 A796. 



£ (xi£i,x 2 £ 2 ) =i J dx 3 dx 4 dt 3 di 4 VF (xit^, x 3 f 3 ) G (xiii,x 4 £ 4 ) V (x 4 i 4 , x 2 t 2 , x 3 t 3 ) 
14 7 (xiti, x 2 f 2 ) =u (ri, r 2 ) + J dx 3 dx 4 dt 3 dt 4 VK (xi*i, x 3 i 3 ) P (x 3 t 3 , x 4 t 4 ) u (x 4 t 4 , x 2 t 2 ) 

P (xiii, x 2 i 2 ) = - i J dx 3 dx 4 di 3 dt 4 G (x 2 i 2 , x 3 i 3 ) G (x 4 f 4 , x 2 i 2 ) T (x 3 i 3 , x 4 t 4 , xifi) 

T (xi*i, X 2 *2, x 3 t 3 ) =5 (xi - x 2 ) S (x 3 - x 4 ) 5 (*! - t 2 ) 5 (ti - t 2 ) 

<5£ (x 1 t 2 ,x 2 t 2 ) 



(5.49) 

(5.50) 

(5.51) 
(5.52) 



/ 



+ / dx 4 dx5dx 6 dx7dt 4 d£5dt 6 dt7- 



SG (x 4 i 4 ,x 5 t 5 ) 
x G (x 4 i 4 , x 6 i 6 ) G (x 7 i 7 , x 5 i 5 ) T (x 6 t 6 , x 7 t 7 , x 3 i 3 ) 

With 

• T4 7 (xiti, x 2 t 2 ) - the screend Coulomb interaction 

• P(xiti,x 2 t 2 ) - the Polarisability 

• T (xiti, x 2 i 2 , x 3 i 3 ) - the vertex function 

This is an exact set of equations that expresses the many body problem in 5 quantities G, P, W, £ and 
r. The complexity of the full many-body problem is shifted to the vertex function. 
=>■ find approximations for V 

A simpe approximation is to neglect the 2nd part of T i.e. T := I 
=>■ Hedin's GW approximation 



£(xi,x 2 ,t) = iG(xi,x 2 ,t) W(xi,x 2 ,t + ) 
P(xi,x 2 ,i) = -iG(xi,x 2 ,t)G(x 2 ,xi,-f) 
W stays as it is 



(5.53) 
(5.54) 
(5.55) 



The real merit of Hedin's equations lies in the fact that the bare Coulomb interaction has been replaced by 
the screened Coulomb interaction. This is important for solids and systems with large polarizability. 

W (r, r', t)= j dr'^'e- 1 (r, r", t-t')v (r" - r) (5.56) 

e (r, r', t) = 5 (r - r') - J dr"dt'v (r - r") P (r", r', t - t') (dielectric function) (5.57) 

The procedure of a typical GW calculation is: 



62 



5 Green 's function theory 

1. Perform a DFT calculation 

2. Build Kohn Sham Green's funtion Go 

3. Calculate Po — —IGqGq 

4. Calculate e = 1 — vP 

5. Invert e 

6. Calculate Wo = e~ 1 v 

7. Calculate a = iG W 

8. Solve Dyson's equation which can be rewritten as 

(Mr) + v H (r)) Mr) + J dr% (r, r', e° w ) ^ s (r') = ef w ^ s (r) (5.58) 



A typical approximation is ^> s (r) = ip^ s (r) 




l r xl r xl r x 

Figure 5.4: Band structure of silicon. 
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Figure 5.5: Band gaps of selected compounds compared to experimental values. 
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